Small-sample corrections for score tests in 
Birnbaum-Saunders regressions 
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Abstract 

In this paper we deal with the issue of performing accurate small-sample 
inference in the Birnbaum-Saunders regression model, which can be useful 
for modeling lifetime or reliability data. We derive a Bartlett-type correction 
for the score test and numerically compare the corrected test with the usual 
score test, the likelihood ratio test and its Bartlett-corrected version. Our 
simulation results suggest that the corrected test we propose is more reliable 
than the other tests. 

Key Words: Bartlett-type correction, Birnbaum-Saunders distribution, chi- 
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1 Introduction 

The Birnbaum-Sau nders distribution, also known as the fatigue life distribution, 
was introduced by iBirnbaum and Saunders! (119691 ) and has received considerable 
attention in recent years. It was originally derived from a model for a physical 
fatigue process w here dominant crack growth causes failure. It was later derived by 
Desmondl (119851 ) using a biologi cal model which followed from r elaxing some of the 
assumptions originally made by lBirnbaum and Saunders! (19_6_9j ). 

The random variable T is said to have a Birnbaum-Saunders distribution with 
parameters a,i] > 0, say B-S(a,r)), if its cumulative distribution function (cdf) is 
given by 

a ( v5 ~ 

where $(•) is the standard normal distribution function and a and 77 are shape 
and scale parameters, respectively. This distribution has a number of interesting 
properties: (i) r\ is its median; (ii) for any k > 0, kT ~ B-S(a,kr))] (iii) its hazard 
function equals zero at t = 0, incre ases up to a maximu m value and then decreases 



F T (t) = $ 



t > 0, 



Kundu et all l2008h : (iv) Y = log(T) is sinh- 



towards a given positive level (see 
normal distributed, say Y ~ SJ\f(at, fi, a), with shape, location and scale parameters 
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given by a, = log(rj) and a = 2, respectively (see Section [2]). iRiectd ( 119991 ) derived 
the moment generating function of the sinh-normal distribution and showed that 
it can be used to obtain both integer and fractional moments for the B-S(a,r)) 
distribution. 



Rieck! (119951 ) derived estimators for the parameters of the B-S(a,T]) distribu- 
tion in complete samples and type II symmetrically censored samples. Some in- 
teresting results ab out improved s t atistic al inference for the B-S(a,rj) distribu- 
tion are available in Lemonte et al. (120071). Extensions of the Birnbaum-Saunders 



distr ibution are p r esent ed in ISanhueza et al.l (1 2008) 



Dfaz-Garcfa and Leival (120051 ) 



and iGomes et al.l (120091 ) . iRieck and Nedelmanl (Il99ll ) proposed a log-linear regres- 
sion model based on the Birnbaum-Saunders distribution. Diagnostic t o ols fo r 



the Birnbaum-Sa under s regression model were developed by I Galea et al.l (120041 ) 



Leiva et al.l ( 2007 ) and Ki and Weil (120071 ) . and Bayesian inference was developed 
bv lTisionasI feoOlf ). 

Inference in Birnbaum-Saunders regressions relies on large-sample theory. For 
instance, a x 2 approximation to the distribution of the likelihood rat io statistic can 



be use d to perform approximate inference if the sample size is large. iLemonte et al. 



(120101 ) presented simulation results showing that the (asymptotic) likelihood ratio 
test can be markedly oversized in small or moderate-sized samples. They derived 
Bartlett-corrected likelihood ratio tests that perform much better than the original 
test. 

Rao's score test is an alternative to the likelihood ratio test. It is simpler to use 
since it only requires estimation under the null hypothesis. Also, it i s often less size 
distorted than the likelihood ratio test; see the simulation results in Cordeiro et al. 



(120031 ). Like the likelihood ratio test, the score test uses an asymptotic approxima- 
tion to the null distribution of the test statistic. This approxim ation can be improved 
by ap plying a Bartlett-type correction to the score statistic; see I Cordeiro and Ferrari 
(Il99ll ) for details. The correction needs to be tailored for each application of interest. 
In this paper, our focus is the Birnbaum-Saunders regression model. 

The paper is organized as follows. In Section [2] the Birnbaum-Saunders regres- 
sion model is presented. In Section [3] we present some cumulants of log-likelihood 
derivatives, which are needed in the subsequent section. In Section H] we derive 
a Bartlett-type correction for the score statistic. Numerical results are presented 
and discussed in Section We compare the performance of the test that uses the 
corrected score statistic with tests that use one of the following statistics: score, 
likelihood ratio and its Bartlett-corrected version. Our Monte Carlo simulations 
suggest that the test that uses the corrected score statistic is much more accurate 
than the rival tests in small and moderate-sized samples. Finally, Section closes 
the paper with some conclusions. 



2 The Birnbaum-Saunders regression model 

Let T ~ B-S(a,T]). The density function of Y = log(T) is 



Mv) 



cosh 



aoV27r 



( v9 exp { 



cr \ a 



yeJR 
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and, as before, we write Y ~ SAf(a, fa a). iRieck and Nedelmanl (Il99lh proposed 
the following regression model: 



y i = x i (3 + e i , i = l,2,...,n, (1) 

where yi is the logarithm of the ith observed lifetime, xj = Xj 2 , . . . , Xi P ) contains 
the ith observation on the p covariates (p < n), (3 — fa, ■ ■ ■ , (3 P ) T is a vector of 
unknown regression parameters, and £j ~ SJ\f(a, 0, 2). The distribution of the errors 
£j has the following properties: (i) it is symmetric around zero; (ii) it is unimodal 
for a < 2 and bimodal for a > 2; (hi) its variance is a function of a only, and has 



no closed- form expression, but IRieck! (119891 ) obtained asymptotic approximations for 



both small and large values of a; (iv) oT x Ei converges in distribution to the standard 
normal distribution when a — > 0. 

The log-likelihood function for the parameter vector 6 = (f3 T ,a) T , apart from 
an unimportant constant, can be expressed as 

n n 

*(0)=5>8(6i)-2E& ( 2 ) 

i=i i=i 

where &i = £a(0) = 2a' 1 cosh([y, - fa]/ 2), £ i2 = = 2a" 1 sinh([y, - fa]/ 2) and 

/ij = xj f3, for z = 1, 2, . . . , n. The nxp matrix X = (x\, x 2 , . . . , a? n ) T is assumed 
to be of full rank, i.e., rank(X) = p. 



3 Cumulants of log-likelihood derivatives 

Derivatives of £{6) with respect to the components of f3 and a are denoted by: U r = 
d£(d)/d/3 r , U a = d£(0)/da, U rs = d 2 £(6)/d(3 r d/3 S} U ra = d 2 £(0)/d/3 r da, U rst = 
d 3 £{0)/d[3 r dM(3t, U rsa = d 3 £(0)/d(3 r df3 s da, etc. Further, we use the following 
notation for joint cumulants of log-likelihood derivatives: K rs = K(U rs ), K r>a = 

etc. All k's are assumed to be of order 0{n). By differentiating (j2j) we have 

Ur = \ ^2 Xir ( ^!^ 2 ~ T~ I aIld Ua = + ~ E 

1 i=i ^ yi/ a a . =1 

The score function for f3 is [7/3 = (l/2)X T s, where s = s(0) is an n-vector whose 
ith element £ a £j 2 - W&i- 

As will be seen below, some cumulants of log-likelihood derivatives depend on 

a , Q ={l-erf(^)}exp(A 
H ere, erf (•) represents the e rror f unction. Details on the e rror function can be found 



111 



Gradshteyn and Ryzhikl (120071 ) . For small values of a (lAbramowitz and Stegunl . 
197(1 p. 298) 

a f l-^ + ^l. (3) 



2tt 4 16 
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For numerical evaluation we recommend the use of fl3]) when a < 0.5. 

The complete list of the cumulants of log-likelihood derivatives that are needed 
in the derivation of Bartlett-corrected score statistic in the Birnbaum-Saunders 



regres sion model is given below. Some cumulants were obtained by iLemonte et al. 



( 120101 ) . Using Bartlett identities, we arrive, after long algebra, at K rst = K rSj t 

K r,s,t = 0? 



■ / j -Eir-Eis ; ^r,s,tu ^ / Xi r Xj s XnXj u 

3si, 



and 



'">/-. .s././l ~ 2 I!- 1 ' I S- 1 ' ll 1 ill • 

8=1 



where a liQ , = 2 + 4/a 2 - a , a \/27r/a, s 1)Q , = -2a 2 , a + (s , a - af )Q! )/8, a 2jQ = -{2 + 
7/a 2 - a , aV ^F(l/(2a) + 6/a 3 ) /v 7 ^}^ and s , a = 12 + 2/a 2 + 16/a 4 + a o , Q 0F(l/a + 

\ C 2i j Oi 1 j \f^2i. A.lsO ; ^tolgl ~~ cxct ~~ 



2n _ (2 + a 2 ) A 6n lOn 

a z a 13 z — ' a' 3 a 13 

i=l 

. [ 2(2 + a 2 ) 1^ _/ (2 + a 2 ) )^ 

K. J \ ) i = i 

„3 _ .Aw/1/^,2U1/J 



with a 3ja = 3/a 3 — a , a v 27r(l/(4a ) + 1/a ). Additionally, 



8n 3 [ 2 + a 2 s 2 , a — ai,, 

_ 48? 



a d a z cr 4 



48n f 2 + 11a 2 3(s 2 , a - a 1>a 



n 



1=1 



4 > 1 a 4 1 4o; 2 



, °^*2,a — u>l,a) ST^ 

+ — a wz. 



Xi r Xi S) 



i=l 



where s 2 Q = 6 + 8ao Q "\A"/ (a^V 7 ^) and 04 a = — 10/a 4 — 4/a 6 + ao Q A/7r(a 4 + 10a 2 + 
8)/(a 7 v/2). 



4 Improved Score Tests 

The hypothesis of interest is Hq: (3\ = (3f\ which will be tested against the alter- 
native hypothesis H\: f3 x ^ f3f\ where (3 is partitioned as (3 = (f3j,f3j) T , with 
f3 1 = (/?!, /3 2 , ... , f3 q ) T and (3 2 = (Pq+i, q +2, • • • , P P ) T for q < p. Here, j3f> is a fixed 
column vector of dimension q. When q = p, the null hypothesis is (3 = (3^. 
The score statistic can be written as 

S R = a&a T X 1 (R r R)- 1 Xj8, 

where a\ >a and s were defined above, R = X\ — X 2 (X 2 T X 2 )^ 1 Xj Xi, with the 
matrix X partitioned as X = [X\ X 2 ] when q < p, following the partition of (3. 
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When q = p, X = X\. Finally, tilde indicates evaluation at the restricted maximum 
likelihood estimate under the null hypothesis. 

Our aim is to obtain simple formulas for A±, A 2 and A 3 , the coefficients of order 
rT l that define the Edgeworth expansion of Sr under the null hypothesis. These 
quantities also define the Bartlett-type correction for the score statistic given by 
Cordeiro and Ferraril (119911 ). The Bartlett-type corrected score statistic is given by 



^ k=i ' 



(4) 



with ^ = {A 1 -A 2 + A 3 )/12q, $ 2 = {A 2 -2A 3 )/{12q{q + 2) \ and ft, = A 3 /{12q{q + 
2)(q + 4)}. General expressions for the A's are given by lHarris (Il985h . The co- 
efficients Ai, A 2 and A3 are functions of joint cumulants of log-likelihood deriva- 
tives up to the fourth order. Whenever they depend on unknown parameters, 
they can be evaluated at the restricted maximum likelihood estimates under the 
null hypothesis. The null distribution of is chi-square with approximation er- 
ror reduced from order 0(n~ l ) to 0(n~ 3 / 2 ). The improved statistic S R can be 
nonmonotonic (in the unmodified statistic). Alternative forms of the Bartlett-type 
corrected s core statisti c that are monotonic transformations of Sr can be found in 
Kakizawa (Il996h and I Cordeiro et al. I fll998h . For a detailed survey of Bartlett and 
Bartlett-type corrections in econ ometrics and statistics, the reader is referred to 
Cribari-Neto and Cordeird (119961 ). 

In the Birnbaum-Saunders regression model, a and (3 are globally orthogonal, 
i.e. K rtCt = 0, for r = 1, 2, . . . ,p. This fact yields some simplification in the derivation 
of the A's. In this case, we can write A\ = A^p + Ai t p a , A 2 = A 2j p + A 2t p a and 
A 3 = A 3) p+A 3) p a , where A 1; p, A 2i p and A 3t p are the corresponding A's obtained as if 
a were known; the additional terms, A lt p a , A 2t p a and A 3i p a , represent the contribu- 
tions due to the fact that a is unknown. The additive structure of the A's also holds 



for o ther models, such as the generalized linear models (ICribari-Neto and Ferraril . 
19951 ) and exponential family nonlinear models (IFerrari et all 119971 ). It can be 
shown that A 3t p a = 0. Also, general expressions for A^p a and A 2) p a are given in 
the Appendix. 

Let Z = X(X T X)~ 1 X T , and for q < p, Z 2 = X 2 (X 2 T X 2 )" x X 2 T . Consider 
= Z Z, zf^ = Z d Z d , etc., denoting the Hadamard (elementwise) 
product of matrices and the subscript d indicates that the off-diagonal elements of 
the matrix were set equal to zero. After long algebra, we obtain 



A h p = ffi,«tr{(Z - Z 2 ) d Z 2d }, A 2) p = g 2 , a tr{(Z - Z 2 )f }, A 



3,/3 



0. 



(5) 



where tr(-) represents the trace operator, gi )Ct = — 96si ta /af a and g 2>a — 
When q = we have Z 2 = and, hence, A^p = A 3 ^p = and A 2t p can be written as 
^2,/3 = 92,a^{Z d 2 ^}. Notice that the formulas given in (jSJ) are functions of (Z — Z 2 ) d , 
Z 2d and a; they do not depend on the unknown parameter vector f3. Also, notice 
that A 3 _p = 0. The formulas in ([Sj) involve si mple matrix algebra and can be e asily 
i mplemented in existing software su c h as R (IR Development Core Teaml . 120081 ). Ox 
flDoornikl . l2006h . MAPLE flAbell et all l2002f ). among others. 
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We obtain 



A 12q St \ X A + %3,a , R , 
Ai,p a = UP- q)9A,a + 9B,a + 96,at, A 2 ,/3 a = , (o) 

n n 

where g 3>a = -2\o? s\J a\ a , g> 4jQ = 4(2 + a 2 )s 3ia / (a 2 a 2 la ), g^ a = as 3>a /a 1>a , 
9e,a = -a 2 S4, a /ai, a , s 3jQ = 2(2 + a 2 )/a 3 - a 3j «, s 4 , a = -{4(1 - 2a 2 )/a 4 + a 4 , a }. 
It is interesting to note that A\.p a and given in ([6]) depend neither on the 

model matrix X (except through its rank) nor on the unknown parameter vector 
j3. It is noteworthy that these formulas are very simple and of easy computational 
implementation. 

We now focus on testing the null hypothesis Hq. a = against the alternative 
hypothesis %i. a ^ a^°\ where is a known positive constant. The score statistic 
for testing can be written as 

& = £(6-i) a , 

where £2 = T^^i^/n. After some algebra, we have 
24n 

Ax = 4 F 2 {(2 + a 2 ) 2 (p + 6) -Aa 3 (2 + a 2 )a 3 , a -a 2 (4 + 5a 2 )a ha } 1 



12 f A(2 + a 2 )p] , „ 40 

A 2 = —{3- ! -\ and A 3 = — . 

n y a di, a J n 

Again, the formulas for the A's are very simple, depend on X only through its 
rank and do not depend on the unknown parameter (3. Clearly, the A's should be 
evaluated at a^°\ 



5 Monte Carlo simulations 

We shall now report Monte Carlo simulation results on the finite sample behavior of 
four tests in the Birnbaum-Saunders regression model, namely: the likelihood ratio 
test (LR), the Bartlett-corrected likelihood ratio test (L-Rb)Jl]the original score test 
(Sr) and the Bartlett-corrected score test We consider the model 

Hi = Pixa + (3 2 Xi 2 H h (3 p x ip + Si, 

where xn = 1 and £j ~ SM(a, 0, 2), % = 1,2, ...,n. The covariate values were 
selected as random draws from the U(0, 1) distribution. The number of Monte 
Carlo replications was 10,000, the nominal levels of the tests were 7 = 10% and 
5%, and all sim ulations were carried out using the Ox matrix programming language 



flDoornikl . 120061 ) 



At the outset, the null hypothesis is = (3 p = 0, which is tested against 

a two-sided alternative, the sample size is n = 25 and a = 0.5 and 1.0. Different 
values of p were considered. The values of the response were generated using f3\ = 
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Table 1: Null rejection rates (%); a = 0.5 and 1.0, with n = 25. 



a = 0.5 







7 = 


10% 






7 = 


5% 




p 


LR 


LR b 


Sr 


c* 
D R 


LR 


LR b 


Sr 


°R 


3 


12.96 


9.99 


9.89 


9.55 


6.96 


5.16 


4.66 


4.91 


4 


14.26 


10.71 


11.35 


10.17 


7.80 


5.18 


5.17 


4.81 


5 


15.63 


10.75 


12.57 


9.93 


8.79 


5.39 


5.89 


4.87 


G 


17.04 


11.15 


13.74 


10.24 


10.10 


5.99 


7.12 


5.11 


7 


19.02 


11.94 


15.38 


10.54 


11.40 


6.26 


7.83 


5.38 


8 


20.73 


12.54 


17.21 


10.83 


12.93 


6.81 


9.54 


5.43 


9 


22.76 


13.30 


18.77 


10.84 


14.41 


6.88 


10.59 


5.22 


10 


24.83 


14.42 


21.19 


11.40 


16.46 


8.20 


12.41 


5.98 










a = 


1.0 












7 = 


10% 






7 = 


5% 




P 


LR 


LR b 


Sr 


°R 


LR 


LR b 


Sr 


a* 
°R 


3 


12.82 


10.04 


9.00 


9.77 


6.79 


5.18 


4.11 


4.91 


4 


14.00 


10.68 


9.88 


10.02 


7.67 


5.20 


4.27 


4.77 


5 


15.41 


10.91 


10.92 


10.17 


8.49 


5.59 


5.06 


5.17 


G 


16.76 


11.29 


12.27 


10.53 


9.76 


5.96 


5.82 


5.23 


7 


18.41 


11.93 


13.52 


11.08 


10.93 


6.29 


6.54 


5.51 


8 


20.65 


12.54 


15.60 


10.92 


12.68 


6.79 


8.01 


5.65 


9 


22.18 


13.43 


17.00 


11.35 


14.36 


6.86 


8.90 


5.41 


10 


24.65 


14.61 


19.43 


11.85 


16.09 


8.22 


10.93 


6.23 



02 = • • • = ftp-2 = 1- The null rejection rates of the four tests are presented in 
Table ffl 

Note that the likelihood ratio (LR) and the score (Sr) tests are markedly liberal; 
the size distortion increases with the number of regressors. For instance, when p = 7, 
a = 0.5 and 7 = 10%, their rejection rates are 19.02% and 15.38%, respectively. It 
is noticeable that the score test is less liberal than the likelihood ratio test. Both 
corrected tests (LRb and S£) are much less size distorted than the uncorrected 
tests. The best performing test is the Bartlett-corrected score test (S£); it displayed 
rejection rates closer to the nominal levels in all cases. In the situation mentioned 
above, the null rejection rates of the corrected tests are 11.94% (LRb) and 10.54% 
(S£). Similar results hold for a = 1.0. More importantly, simulations carried out for 
a wide range of values of a reveal an analogous pattern (results not shown here for 
the sake of space) H For instance, when n = 25, p = 4 and 7 = 10%, the rejection 
rates of H : /3 3 = /3 4 = are 14.33% (LR), 11.42% (LR b ), 11.76% (S R ) and 10.08% 
(££) for a = 0.1, and 12.44% (LR), 11.23% (LR b ), 5.12% (S R ) and 9.52% (Sg) for 
a = 10. We noticed that the uncorrected score test becomes conservative when a is 
large. 

Table |2] reports results for sample sizes ranging from 15 to 100, a = 0.5 and 
p = 7. The null hypothesis under test is Ho'- = Pi = 0. The uncorrected tests are 



1 LR h = LR/{1 + B(0)/q}, where B(6) is the Bartlett factor; see iLemonte et al.l (|2010l) for 
details. 

2 We used the approximation in ([3]) for ao jQ when a < 0.5. 
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remarkably oversized when the sample size is small. As expected, the null rejection 
rates of all tests approach the corresponding nominal levels as the sample size grows. 
Clearly, the corrected tests are less size distorted than the unmodified tests, and the 
test that uses S R has superior behavior than the tests that use LR, LR b and Sr in 
all cases. For example, when n = 20 and 7 = 5%, the null rejection rates are 13.87% 
(LR), 7.06% (LRb), 9.15% (S R ) and 5.62% (£*). 

Table 2: Null rejection rates (%); a = 0.5, p = 7 and different sample sizes. 



n 




7 = 


10% 






7 = 


5% 




LR 


LRb 


Sr 


c* 
°R 


LR 


LRb 


Sr 


c* 


15 


27.93 


16.03 


21.27 


12.96 


19.22 


9.44 


12.01 


7.32 


20 


21.93 


13.03 


17.29 


11.10 


13.87 


7.06 


9.15 


5.62 


30 


17.28 


11.55 


14.53 


10.56 


10.13 


5.69 


7.51 


5.16 


40 


15.26 


10.87 


13.19 


10.45 


8.59 


5.58 


6.78 


5.18 


50 


14.10 


11.05 


12.74 


10.75 


8.18 


5.59 


6.64 


5.34 


100 


11.55 


9.97 


10.81 


9.87 


6.03 


5.01 


5.48 


5.00 



We shall now present simulations results on the powers of the tests. We set 
p = 3, a = 0.5 and n = 30 and 50. Since the unmodified tests are considerably 
oversized, we have only considered the two corrected tests. The rejection rates were 
obtained under the alternative hypothesis 02 = 03 = 8, with different values of S 
(5 > 0). Note that the two tests display similar powers. For instance, when n = 50, 
7 = 10% and 5 = 0.5, the nonnull rejection rates are 86.90% (LR b ) and 86.70% 
(S£). Not surprisingly, the powers of the tests increase with n and also with 5; see 
Table EJ 

Table 3: Nonnull rejection rates (%); a = 0.5, p = 3 and n = 30 and 50. 



LR b 


Sr 


n S 10% 5% 


10% 5% 


30 0.1 12.69 6.81 
0.3 31.68 20.90 
0.5 64.43 50.32 
0.7 88.54 80.50 


12.47 6.62 
31.29 20.30 
63.64 49.26 
88.26 79.68 


50 0.1 14.80 8.03 
0.3 48.87 35.28 
0.5 86.90 78.48 
0.7 99.03 97.54 


14.76 7.97 
48.65 34.88 
86.70 78.17 
99.05 97.50 



We performed Monte Carlo simulations considering hypothesis testing on a. To 
save space, the results are not shown. We noticed that the corrected tests display 
superior behaviour than the uncorrected tests. For example, when n = 30, p = 4, 
7 = 10% and "H : a — 1-0, we obtained the following null rejection rates: 19.86% 
(LR), 10.73% (LR b ), 13.77% (S R ) and 9.89% Again, the best performing test 

is the Bartlett-corrected score test. 
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As noted by the referee, our results can also be used to inv ert the Edgeworth 



expansion of the sco r e stat istic null distribution function; see lHarrisI (jl985[ ) and 



Cordeiro and Ferraril (119911 ) for details. The test is then performed by comparing 
the unmodified statistic with a corrected critical value. To be specific, let 7 be the 
desired level of the test and gi_ 7 be the 1 — 7 quantile of the x 2 limiting distribution 
of the test statistic. The corrected critical value is 

9l % = gi- 7 (l + X;^^), (7) 
^ k=i ' 

where for i = 1, 2, 3, are given just below (j3J). We shall now present some simula- 
tion evidence on the finite sample behavior of the score test that uses the corrected 
critical value {Sh)- For the sake of comparison, the new set of simulation results 
includes rejection rates of the other four tests {LR, LR b , Sr and S R ) and the score 
test that uses parametric bootstrap critical values {S 000 t)- The bootstrap correct test 
is performed as follows. First, one generates B bootstrap samples (we set B = 600) 
from the assumed model with the parameters replaced by restricted estimates com- 
puted using the original sample, under T-L , i.e. imposing the restrictions stated in 
the null hypothesis. Second, for each pseudo- sample, one computes the score statis- 
tic; S R denotes the score statistic for the 6-th sample, b = 1,2, . . . , B. Third, the 
1 — 7 percentile of Sr is estimated by <fi_ 7 , such that #{S R < <fi- 7 } / B = 1 — 7. 
Finally, one rejects the null hypothesis if Sr > §i_ 7 . 

Table H] reports rejection rates of the different tests of the null hypothesis %q\ 
f3 s = f3 9 = for a = 1.5, n = 30 and p = 9. The score test that uses the corrected 
critical value ([7]) and the Bartlett-corrected score test have similar size properties and 
are less size distorted than both likelihood ratio tests (corrected and uncorrected) 
and the original score test. The test that uses the bootstrap corrected critical value 
performs very well although being slightly conservative. Its main disadvantage over 
the other tests is the need of a computer intensive procedure. 



Table 4: Null rejection rates (%); a = 1.5, n = 30, p = 9. 





7 = 10% 


7 = 5% 


LR 


19.36 


11.71 


LRb 


12.17 


6.40 


S R 


14.79 


7.62 


S£ 


11.25 


5.66 


Sh 


11.42 


5.64 


Shoot 


9.53 


4.71 



Finally, we end this section by reporting the first four moments of LR, LR 0} Sr 
and S R and the corresponding moments of the limiting \ 2 distribution in the setting 
of Table HI see Table [5j Clearly, the best agreement between the true moments 
(obtained by simulation) and the moments of the limiting distribution is achieved 
by the Bartlett-corrected score statistic, S* R . 
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Table 5: Moments; a — 1.5, n — 30, p — 9. 

LR LRi, Sr Sft Xi 
Mean 2^81 2^20 ~2A\ 2^10 2.00 

Variance 8.08 4.94 4.89 4.32 4.00 
Skewness 2.17 2.17 1.59 1.98 2.00 
Kurtosis 10.93 10.88 6.46 9.19 9.00 



6 Conclusions 

We addressed the issue of performing testing inference in Birnbaum-Saunders re- 
gressions when the sample size is small or moderate. We derived a Bartlett-type 
correction for the score test. We numerically compared the behaviour of four tests, 
namely: the likelihood ratio test (LR), its Bartlett-corrected version (LR b ), the 
score test (S R ) and its Bartlett-type corrected counterpart (<S£). Also, we perfomed 
a small simulation study that includes two score tests that use critical values that 
are approximated either analytically or by a parametric bootstrap method. The 
bootstrap-corrected test performs well, but requires a computer intensive procedure. 
Overall, our numerical results favor the test obtained from applying a Bartlett-type 
correction to the score test statistic. Therefore, we recommend the corrected score 
test proposed in the present paper for practical applications. 
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Appendix 



Since K r ,a = 0, for r = 1,2, ... ,p, the total Fisher information matrix for 6, K(G), and 
its inverse, K(6) _1 , are block diagonal. Using the partition of X induced by the null 
hypothesis 7-Lq: Pi = (3^ , we can write 





~K U 


K 12 


" 




~K n 


R 12 





K(6) = 




K 22 





and K(6)- 1 = 




K 22 


























Let 




K 



22 



K 



R 21 



K 12 

R 22 



A3 





and M 



M/3 




where M/3 = — Ap. Let m ra and a ra be the (r,p + l)-th elements of M and 

A, respectively, and let m aa and a aa be the (p + l,p + l)-th elements of M and A, 



respectively. We have m r 



0, a r 



(r = 1,2, ... ,p) and 
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a aa = K> a a . Therefore, the general formulas for A\^ a and Aifi a are 

A\fi a = 3 y {l^aak 4 2rCo,afc)( K rst + lK rs ^t)a a a a st' m kr 

El 
\K'ijk 4 < ^K'i,jk)\K'rao( 4 ^^ , ra,a)^'ij^'aa^'kr 

4 3 y ' (ftoaJ; "I - 2Avq, iQ ,^) (/tract "I - 2K J - Q . jCe )ci Q , Q ,T7l^ r 

6 ^ ^ {^ija 4 ^^ , i,ja)^'a,s,t^'ij^'aa^'st ^ ^ ~] (^aaa 4 2'C Q;jQ , Q ,)K a!jSj j CL aa 7Jl s i 
4 6 / { K i,ak ~ K i,a,fc)( K rat 4" '^l^ra,t)0'aa a 'kt' m ir 

El 

4 6 ^ ^ (^i,C(ce ( ^rcea 4 2K r .o,a)^oa:^''i'' 6 ^ ^ {^i,j,a,a 4 l^iJ.aa)0'aa^T'ij 

and 

^2,/3a = ~3 ^ ^ Ki.j,a K a,s,t a aa m ij m st 4 6 ^ ^ ( K a«fc 4 2K aa i^K rs j(laa' m kr rn 'st 

where ^ denotes summation over all the components of /3, and ajj e m,,,- = 1,2, ... ,p) 
are the (z,j)-th elements of A and M, respectively. 

Plugging the cumulants given in Section 0] in the formulas for A\^ a and A2^ a , we ob- 
tain ([6]). To save space, we will only show how to obtain —6 Y' ( K ija+2Ki.j a )K ajS: taija aa m s t. 
The other terms can be obtained in a similar fashion. Notice that 

aX^' ( , 9 \ _ 6(2 + a 2 )s 3ta a aa \ ^ 

J y^ija 4" ^l^i,ja)^a,s,tO'ijO'aa1^T'st — 3 / / XiiXljX ms X m fCLij7Jl s t. 

l,m=l 

Inverting the order of summation and rearranging the terms, we have 
-rX^'( x9 ^ _ 3(2 + a 2 )g 3i o A w w 

l,m=l 

The terms Yl' a ij x n x ij an d Yl m stx ms x m t represent the elements (/,/) and (m,m) of 
4^2/01,0 and 4(Z — Z<i)la\ t0l , respectively. The matrices Z = {z/ m } and Z2 = {z2i m } 
were defined in Section [H Hence 



R W. n _ 48(2 + q 2 )^ A , 

O / \Kija ~r ^ K i,jaj K a,s,t a ij a aa m st — 2 / ^ z 2U\ z mm z 2r 



naa\ 

x '" (,m=l 



Rom S"m=l z 2u{zmm - z 2m m) = tr(Z 2 )tr(Z - Z 2 ) = (p - g)<?, we obtain 
-fiST'f x9 ^ _ 48g 3 , Q (2 + a 2 )(p - q)q 
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